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ABSTRACT 

We report on the first global three-dimensional (3D) MHD simulations of disk accretion onto a rotating 
magnetized star through the Rayleigh- Taylor instability. The star has a dipole field misaligned relative 
to the rotation axis by a small angle Q. Simulations show that, depending on the accretion rate, a star 
may be in the stable or unstable regime of accretion. In the unstable regime, matter penetrates deep 
into the magnetosphere through several elongated "tongues" which deposit matter at random places 
on the surface of the star, leading to stochastic light-curves. In the stable regime, matter accretes in 
ordered funnel streams and the light-curves are almost periodic. A star may switch between these two 
regimes depending on the accretion rate and may thus show alternate episodes of ordered pulsations and 
stochastic light-curves. In the intermediate regime, both stochastic and ordered pulsations are observed. 
For > 30°, the instability is suppressed and stable accretion through funnel streams dominates. 



1. INTRODUCTION 

Magnetospheric accretion occurs in different stars, for 
example in classical T Tauri stars (CTTSs), which are the 
progenitors of Solar-type stars (e.g., Hartmann 1998; Bou- 
vier et al. 2007), in magnetized white dwarfs in some cat- 
aclysmic variables (e.g., Warner 1995; Warner & Woudt 
2002), and in accreting millisecond pulsars, which are 
weakly magnetized neutron stars (e.g., van der Klis 2000). 
The matter close to the star may either flow around the 
magnetosphere, forming ordered funnel streams going to- 
wards the magnetic poles of the star (e.g., Ghosh & Lamb 
1978; Camenzind 1990; Konigl 1991), or may accrete di- 
rectly through the magnetosphere due to the Rayleigh- 
Taylor (RT) instability (e.g., Arons & Lea 1976; Eisner & 
Lamb 1977; Scharlemann 1978). Earlier 2D and 3D nu- 
merical simulations have shown accretion through funnel 
streams (Koldobaet al. 2002; Romanova et al. 2002, 2003, 
2004). Recent 3D MHD simulations performed for a wider 
range of parameters have shown that in many cases the 
disk-magnetosphere boundary is RT unstable. 

Previous theoretical models and non-global simulations 
gave useful but restricted analyses of this problem. Arons 
and Lea (1976) investigated magnetospheric accretion 
through the RT instability in the non-rotating case with a 
spherical accretion geometry. Spruit and Taam (1990) in- 
vestigated the stability of an infinitely thin rigidly rotating 
disk. The RT instability in magnetized disks was studied 
by Kaisig et al. (1992); Spruit et al. (1995); Lubow & 
Spruit (1995). Li and Narayan (2004) investigated disk- 
magnetosphere interaction in the case of an infinitely thick 
disk with a vertical magnetic field. 2D simulations have 
been performed by Wang and Nepveu (1983) and Wang 
& Robertson (1985), while Rastatter and Schindler (1999) 



performed both 2D and 3D simulations, but in a patch 
(see also Stone and Gardiner 2007). Such simulations, 
while shedding light on many important features of the 
instability, do not take into account many factors which 
are present in global simulations, such as the possibility 
of matter fiowing through funnel streams to the magnetic 
poles of the star. 

In this paper we show results from global 3D MHD simu- 
lations, where the simulation region includes the disk and 
the whole magnetosphere of the star. This paper sum- 
marizes these new results, while in the following paper 
(Kulkarni & Romanova 2008, hereafter - KR08) we present 
results for a larger range of parameters. 

2. MODEL AND RESULTS OF 3D MHD SIMULATIONS 

2.1. Model. We use a "cubed sphere" Godunov-type 
numerical code (Koldoba et al. 2002) and solve the full set 
of 3D MHD equations with initial and boundary conditions 
similar to those used in Romanova et al. (2004). Simu- 
lations were done for non-relativistic as well as relativis- 
tic neutron stars. We approximate relativistic effects us- 
ing the Paczynski-Wiita potential (Kulkarni & Romanova 
2005). Our model includes viscosity to regulate the matter 
flux M through the disk. The viscous stress in the disk is 
proportional to the a-parameter, which we varied in the 
range a = 0.02 — 0.3. Most of the results shown here are 
for the grid Nr x x Ny — 148 x 61 x 61 in each of the 
6 blocks of the "cubed sphere" . Test simulations for twice 
as fine a grid and twice as coarse a grid show that the 
number of modes does not depend on the grid resolution. 

The simulations are done in dimensionless form and are 
applicable to stars over a wide range of scales, if the mag- 
netospheric radius is not very large compared to the 
radius of the star R*, r^ = (4 — 5)i?*; is determined by 




Fig. 1. — (a) An example of stable accretion. The surface is a constant density surface, and the lines are sample magnetic 
field lines. The magnetic axis is misaligned relative to the rotational axis hy Q = 15°. (b) An example of unstable 
accretion for = 5°. The lower panels show light curves from the hot spots. Time is measured in orbital periods at 
r = 1. Only the inner part of the simulation region is shown. 



the balance between the magnetospheric and matter pres- 
sure, so that the modified plasma parameter at the disk- 
magnetosphere boundary /3 — {p+ pv^)/ {B^ /8n) « 1. The 
reference units are as follows: the length Rq = R^,/0.35, 

velocity Vo — (GM*/i?o)2, period Pq = 27ri?o/i'Oj den- 
sity pq = Bq/vq, magnetic moment pq ~ BqRq, accretion 
rate AIq — pqVqRq (see KR08 for a complete description 
of units) . We show results for a star with a dimensionless 
magnetic moment p — 2. We varied the period P» of the 
star so that the corotation radius Vcor — {GM/fliy^-^ in 
dimensionless units varied in the range rcor = 1.2 — 3 (with 
Tcor ~ 2 in the main case). The initial disc structure is 
determined by the initial disk-corona-magnetic field equi- 
librium (e.g., Romanova et al. 2002) and it has the same 
fiducial density in all runs, so that the accretion rate in 
the disk is regulated mainly by the a-parameter of viscos- 
ity. The disk is sufficiently large {Rd ~ 48i?*) to supply 
matter during the whole simulation run, which lasted up 
to 50 periods of rotation Pq. Most of the runs performed 
were for the small misalignment angle = 5° which helps 
excite the perturbations. Runs for larger O have shown 
that the instability is present up to O « 30° (sec KR08). 
2.2. Stable and unstable regimes of accretion. In 
our earlier work we performed 3D MHD simulations of 
the disk-magnetosphere interaction for a variety of mis- 
alignment angles from O = 0° to O = 90°, and showed 
that matter flows from the inner regions of the disk to the 
star in symmetric funnel streams, which are quasi-stable 
features over long times (Romanova et al. 2003, 2004). 
Figure la shows an example of accretion through funnel 
streams for O = 15°. In this type of flow, matter is lifted 
above the equatorial plane and then flows along field lines. 
The str cams hit the star near the magnetic poles, forming 
ordered hot spots. 



Next we chose a case with a small misalignment angle 
(0 = 5°) and increased the accretion rate in the disk, in- 
creasing the a-parameter. We observed that at a critical 
value of a = 0.04, the instability appears and the accretion 
acquires an essentially different pattern: the RT instability 
develops at the inner edge of the disk, and matter accretes 
through equatorial "tongues" which penetrate deep into 
the magnetosphere. The structure of the tongues is oppo- 
site to that of the funnel streams: they are narrow and tall 
(see Figure lb). Such a tongue shape is determined by the 
geometry of the magnetic field and has been predicted by 
Arons and Lea (1976). Simulations performed for different 
parameters and grid resolutions have shown that the num- 
ber of growing modes is always small, m = 2 — 7 (KR08). 
The matter energy-density dominates inside these tongues. 
Test runs at twice as large density in the disk have shown 
that the instability appears even at a = 0.02 which proves 
that the matter accretion rate, and not a, is important. 
During the interchange process, the heavy ffuid elements 
of the disk (matter dominated) change positions with the 
light ffuid elements of the magnetosphere (magnetically- 
dominated). The magnetospheric plasma carries frozen- 
in magnetic field lines, so that displacing the light ffuid 
leads to the effect of pushing magnetic field lines aside 
(see Figure 2). Equatorial shces of the density distribu- 
tion for two moments of time are shown in Figure 3 (see 
also http:/ /astro. cornell.edu/us-rus/stereo. htm for anima- 
tions). 

2.3. Light-curves from the hot spots on the surface of 
the star were calculated, to study the difference in observa- 
tional properties between the stable and unstable regimes. 
We assume that the kinetic energy of the stream is con- 
verted at the surface of the star into isotropic black-body 
radiation. Integration of radiation in the direction of the 
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Fig. 2. — Penetration of the tongues through the magnetosphere. The tongues are shown by density contours in the 
equatorial plane (the highest density is shown in red, and the lowest in dark-blue, with a density contrast of about 300). 
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Fig. 3. — Density distribution in the equatorial plane after 20 and 30 periods of rotation at ?■ = 1. The background shows 
the density (the colors have the same meaning as in Figure 2). The black line shows the position of /? = 1. The plots are 
shown in the reference frame rotating with the star. 



observer is performed (see details in R04). In the stable 
regime, the position of the hot spots is almost constant, 
so that the light-curve is almost sinusoidal (see Figure la, 
bottom panel). In the strongly unstable regime, sporadi- 
cally forming tongues hit the star in random places, so that 
the light curve is irregular (see Figure lb, bottom panel). 
In the intermediate cases both funnels and tongues are 
present, and mixed light-curves are expected. If a certain 
number of tongues dominates (for example, the m = 2 
mode dominates in some cases) then quasi-periodic os- 
cillations are expected. A more detailed analysis of the 
light-curves in different cases will be presented in a future 
paper. 

2.4. Onset of instability. We compared our simula- 
tions with a few relevant theoretical approaches. A heavy 
fluid supported against gravity by a lighter one is RT un- 
stable unless there is some opposing force. A homoge- 
neous vertical field at the disk-magnetosphere boundary 
does not oppose the growth of i/i-modes (see e.g. Chan- 



drasekhar 1961). A more general criterion for differen- 
tially rotating magnetized accretion disks indicates that 
the disk is unstable to excitation of </)-modes if 7^j^ = 
-geffd\n{Y,/Bz)/dr > 2{r,ndQ / dr)'^ = jf^ (Spruit et al. 
1995; Spruit & Taam 1990; Kaisig et al. 1992; Lubow & 
Spruit 1995). Here, —geff = — f^^) is the effective 

gravity, and S = 2 phis the surface density. That is, for the 
instability to start, the surface density per unit magnetic 
field strength T./Bz should drop off fast enough in the di- 
rection of the effective gravity (—geff), so that the term 
7^2 is larger than the term associated with the shear, 7^, 
which tends to oppose the instability by smearing out the 
perturbations. We calculated the (/)-averaged values of 7^^ 
and 7q for the stable and unstable cases (see Figure 4) and 
compared the positions of the curves in the inner regions 
of the disk where the instability starts, that is at r ^ r„j. 

One can see that in the stable case (a — 0.02), 7^^ ^ 7a 
while in the unstable case {a — 0.1) » 7^. Thus the 
theory makes a reasonably good prediction for the onset of 
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Fig. 4.— Parameters corresponding to the stability criterion by Spruit et al. (1995) described in §2.4. 
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Fig. 5. — Radial distribution of the density p, angular velocity and Keplerian angular velocity Q.k in the equatorial 
plane. The dashed line marks the position of the magnetospheric radius where /? = 1. 



instability. The term 7^^^ varies strongly mainly because 
in the unstable case, a significant amount of matter accu- 
mulates near and the density drop-off is sharper than 
in the stable case (see Figure 5). 

The RT modes observed in simulations are excited by 
different mechanisms. One of them is the small non- 
axisymmetry associated with the misalignment of the 
dipole. Another is the inhomogeneities in the inner regions 
of the disk driven by the azimuthal component of the field 
trapped inside the inner regions of the disk. The Kclvin- 
Helmhotz instability seems to be less significant (see also 
Rastatter & Schindler 1999). 

2.5. Why low-m modes dominate. According to the 
theory of the RT instability, high-m modes should grow 
faster than low-m ones (e.g. Chandrasekhar 1961). In a 
rotating disk, the shear dSl/dr may efRciently damp high- 
m modes (e.g., Lubow & Spruit 1995) and may be the 
main reason for damping of the high-m modes in our sim- 
ulations. Initially, in high-resolution simulations we see 
the formation of m = 20 — 30 modes, but in one rotation 
period they are damped, and only the low-m modes grow. 
KR08 analyzed the number of modes based on the Li & 
Narayan (2004) criterion and reached a similar conclusion. 

Another possible reason for the suppression of the high- 
m modes is the presence of the azimuthal component of 
the field at the inner edge of the disk which typically con- 
stitutes ~ (5—30)% of the poloidal component. The damp- 
ing effect of is expected from theory (Chandrasekhar 
1961), and has been noticed by Wang & Nepveu (1983) 
and Rastatter & Schindler (1999). 

2.6. Boundary between the stable and unstable 
regimes. We performed multiple simulation runs for iden- 
tical initial density distributions in the disk but different 



viscosity parameters a and dimcnsionless periods of the 
star p^*™ (these bulk simulations were done for a grid 
72x31^). The symbols in Figure 6a correspond to differ- 
ent runs. One can see that there is a boundary between 
the stable and unstable cases with unstable cases corre- 
sponding to higher a. For similar initial disks and ac- 
cretion rates through the disks, M is proportional to a, 
and thus the instability occurs at high enough accretion 
rates through the disk. The boundary is higher at smaller 
periods P^^"^, because at faster rotation of the star the ef- 
fective gravity becomes less negative and thus less matter 
may accumulate at the inner edge of the disk. That is why 
a larger a is required to bring more matter to the inner 
edge of the disk to satisfy the instability criterion. Cases 
with periods P^*™ w 2.7 — 3 correspond approximately to 
the rotational equilibrium state (Long et al. 2005). Stars 
with smaller periods spin down and, at P_^™ ^ 1.8, enter 
the "propeller" regime, where rcor < i"m (e.g. Illarionov 
& Sunyaev 1975; Lovelace et al. 1999; Rappaport et al. 
2004). In this regime a star tends to expel a significant 
part of incoming matter (Romanova et al. 2005; Ustyu- 
gova et al. 2006). At high enough a we observed events of 
instability even in this regime (see also Wang & Robertson 
1985), but the simulation runs were brief, and this region 
requires further analysis. 

We calculated the accretion rate onto the surface of the 
star M^*™ for the above runs. The boundary between sta- 
ble and unstable accretion M*™ w 0.3 seems to be almost 
independent of P^, , which is probably a result of the combi- 
nation of the lower disk density due to the propeller effect, 
and higher a. This boundary is approximate, with inter- 
mediate cases above and below it in which only marginal 
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Fig. 6. — Position of stable (circles), unstable (triangles) and marginal (squares) cases in the (a) a - P^™ and (b) M^™ - 
P^™ planes. The top horizontal axis shows the corotation radius. The solid lines show an approximate boundary between 
stable (below) and unstable (above) runs. The dash-dotted line marks the approximate boundary of the propeller regime. 



instability is observed. In this transition region, both fun- 
nels and tongues are present. Higher up, for M^™ > 0.5, 
the instability is strong and accretion through tongues 
dominates. We note that at low enough accretion rates 
through the disk even stable accretion may shut-off, be- 
cause there is insufficient pressure gradient at the disk- 
magnetosphere boundary to form magnetospheric funnel 
streams. 

2.7. Application to different magnetized stars. 

The results of our simulations can be applied to a vari- 
ety of magnetized stars with small magnetospheres, ~ 
(4 — 5)i?*, including CTTSs, weakly magnetized white 
dwarfs (WDs) in cataclysmic variables and neutron stars 
(NSs) in LMXBs. The dividing line in Figure 6b corre- 
sponds to M^™ Ri 0.3. Below we show the accretion rate 
corresponding to this boundary in dimensional units for 
different types of stars: 

CTTSs : Mr « 2.1 x 10"**(P*/10^G) Vim-iMo/yr, 



WDs : M," « 1.4 X lO-^(S,/lO^G)Vim-iM0/yr, 



NSs : Mf « 2.2 x 10-^{B^/lO^G)^rim-^ Mq/jt. 
The period in dimensional units is: 

where scaling factors for different stars are: 
CTTSs: Ap « 1.78 days, r = R^/2Rq, m = M*/0.8Mq, 
WDs: Ap « 29.4 sec, r = i?*/5 x 10^ cm, m = M*/1Mq, 
NSs: Ap « 2.22 ms, r = cm, m = M*/1.4M0. 

Note that the dimensional increases with the mag- 
netic field as so that the boundary will be higher for 
higher B^, . This reflects the fact that we have fixed the di- 
mensionless magnetic moment /z, and thus consider stars 
with an approximately fixed ratio r^/i?* « 4 — 5. We 
performed another set of runs for smaller magnetospheres 
(smaller Vm/R* for /i — 0.5) and observed that the di- 
mensionless boundary corresponds to a smaller M*™ (see 



KR08). Note also that all the above runs were done for 
8 = 5°. We expect that for larger the boundary will 
move up, because funnel accretion is more favorable at 
larger G. 

3. OBSERVATIONAL CONSEQUENCES 

The existence of two regimes of accretion changes our 
understanding of accreting magnetized stars and their ob- 
servational properties: the lack of a periodic signal in the 
light curve does not rule out accretion and an ordered mag- 
netic field. For example, the random light-curves in many 
CTTSs or dwarf novae do not preclude an ordered stellar 
field, if the star is in the unstable regime. Depending on 
the accretion rate, episodes of stable and unstable accre- 
tion may alternate, so that periods with pulsations may 
be intermittent and may be followed by periods with no 
pulsations. Recently a number of intermittent accreting 
millisecond pulsars have been discovered (e.g., Kaaret et 
al. 2006; Altamirano et al. 2007; Galloway et al. 2007; 
Gavriil et al. 2007). Although the reason for this behavior 
is not yet understood, we note that in, for example, HETE 
J1900. 1-2455, the pulsations disappeared when the overall 
X-ray fiux increased (Kaaret et al. 2006) which may be an 
example of a transition to the unstable accretion regime. 

In some cases a definite number of tongues may dom- 
inate, which may lead to quasi-periodic oscillations in 
the light curves (Li & Narayan 2004). These oscilla- 
tions may be candidates for one of the QPOs observed 
in Type II (accretion-driven) bursts in LMXBs. Comp- 
tonization of photons by high-energy electrons may lead to 
only a small departure of the light-curve from the thermal 
ones obtained using the approximation of isotropic black- 
body radiation (Poutanen & Gierlihski 2003), so that the 
QPO features may survive Comptonization (see, however, 
Titarchuk et al. 2007). 

In the case of young stars surrounded by gas/dust disks 
where planets are forming and migrating inward, the un- 
stable tongues may support inward migration, so that the 
magnetospheric gap, which can halt migration (e.g. Lin et 
al. 1996; Romanova & Lovelace 2006; Papaloizou 2007), 
may form only in the state of stable accretion. 
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